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SUMMARY 

The conservation-law form of the inviscid gasdynamic equations 
has the remarkable property that the nonlinear flux vectors are 
homogeneous functions of degree one. This property readily permits 
the splitting of flux vectors into subvectors by similarity trans- 
formations so that each subvector has associated with it a speci- 
fied eigenvalue spectrum. As a consequence of flux vector split- 
ting, ne. explicit and implicit dissipative finite-difference 
schemes are developed for first-order hyperbolic systems of equa- 
tions. Appropriate one-sided spatial differences for each split 
flux vector are used throughout the computational field even if the 
flow is locally subsonic. The results of some preliminary numeri- 
cal computations are included. 
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1. INTRODUCTION 


Finite-difference schemes for the conservation-law form of the 
unsteady inviscid gasdynamic equations are restricted to a very 
limited class of spatial difference approximations in subsonic flow 
regions. Only centered difference operators lead to difference 
methods that are simultaneously stable for both the positive and 
negative characteristic speeds (i.e., eigenvalues) that are asso- 
ciated with the spatial flux terms in subsonic flow. Use of any 
other class of spatial differential operator requires splitting 
the flux terms into components of a restricted type. 

There are various reasons for using one-sided spatial differ- 
ence operators. For example, for the model scalar wave equation, 
one-sided (or upwind) schemes frequently have superior dissipation 
and dispersive properties to those of a centered scheme [1,2]. An 
explicit second-order accurate upwind scheme can also have twice 
the stability bound of a centered second-order scheme [1]. Another 
motivation stems from a desire to increase numerical efficiency of 
implicit algorithms. For example, an implicit upwind finite 
difference algorithm can lead to a lower diagonal banded matrix 
that is more easily Inverted than the tridiagonal and pentadlagonal 
matrices usually associated with centered schemes. 

Our objective is to devise a means of splitting the flux 
vectors of a hyperbolic system in order to extend the class of 
allowable spatial differencing schemes to achieve more robust 
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algorithms and to improve computational efficiency. As in earlier 
related work [2,3], we restrict our attention to the inviscid gas- 
dynamic equations in conservation-law form and take advantage of 
the fact that the flux vectors are homogeneous of degree one. We 
have not investigated first-order conservative systems that are 
nonhomogeneous. The basic ideas used here, however, apply to first- 
order nonconservative hyperbolic systems of equations. A related 
explicit algorithm for the equations of gasdynamlcs written in 
nonconservation law form was recently proposed by Moretti [4]. 

In this paper we first review the restrictions placed on the 
spatial difference operators of hyperbolic systems that have both 
positive and negative eigenvalues. Using the one-dimensional 
inviscid equations of gasdynamlcs, we then develop a methodology 
for splitting the equations into components of the same character- 
istic behavior. Both explicit and implicit numerical algorithms 
are devised and tested for the split system of equations. The 
methodology and algorithms are then extended to multidimensions. 

2. MOTIVATION AND BACKGROUND 

In this section we review the restrictions placed on spatial 
difference approximations by the characteristic speeds (eigen- 
values) of a hyperbolic system. 

To illustrate the basic notions we consider a one-dimensional 
system of conservation laws 
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( 2 . 1 ) 


3U . 3F 

+ 0 • 


where U and F are m-component column vectors. The system (2.1) 


can be rewritten as a quasi-linear system 




( 2 . 2 ) 


where A is the Jacobian matrix 3F/3U. The system (2.2) is 
hyperbolic at the point (x,t,U) if there exists a similarity trans- 
formation such that 


Q" *AQ » A * 


(2.3) 


where A is a diagonal matrix, the eigenvalues of A are 

real, and the norms of Q and Q _1 are uniformly bounded. 

For the purpose of a linear stability analysis, we assume 
that the coefficient matrix A is "frozen," that is, constant. 
By virtue of Eq. (2.3), Eq. (2.2) can be transformed to the 
uncoupled system 


3u 3u 

jt ~3x” " ^ ■ 1* 2, 3, . . • , m 


(2.4) 


by defining a new vector u - (uj, u 2 , u 3 , .... u,,,) ■ Q~*U. 

Consequently, when analyzing the stability of numerical algorithms 
as applied to the linearized version of the system (2.2), we need 
only examine the scalar Eq. (2.4). For simplicity, the subscript 
l will be dropped in the remainder of this section. 
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To analyze the effect of one-sided spatial differences on 
stability we leave the tine variable continuous and discretize the 
spatial variable as x ■ ■ jAx. Let 9u/3x be approximated by 

the first-order one-sided difference quotient 


3u 

3x 


j 


Vj 

Ax 


+ 0(Ax) , 


where V x is the classical backward-difference operator 


7 u. 

* j 


“j “ °j-l * 


(2.5) 


( 2 . 6 ) 


This spatial discretization reduces Eq. (2. A) to a system of first- 
order ordinary differential equations: 

u 4 - u 4 


du.i uj 
+ A • 

dt Ax 


j-1 


0 . 


(2.7) 


For simplicity, assume spatially periodic boundary conditions 
and look for a solution of the form 

ikjAx 


Uj (t) » v(t)e 


( 2 . 8 ) 


where v(t) is the Fourier coefficient, i ■ *^T, and k is the 
wave number. By inserting Eq. (2.8) in Eq. (2.7), one finds that 
the Fourier coefficient satisfies the ordinary differential 
equation 


dv 

dF “ 0V ’ 


(2.9a) 


where 


- — ^2 8in 2 0^ + i sin ej , 6 ■ kAx . (2.9b) 
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The solution is v ■ Ce so for Eq. (2.9a) to have a bounded solu- 
tion, the real part of o must satisfy Re a < 0, which requires 
that X > 0. 

Instead of the backward-difference operator (2.6), let 3u/3x 
be approximated by 


Isij ■ h Vj + »«*> - 

where A x is the forward-difference operator 

Vj ■ Vi ' “3 • 

If we repeat the Fourier stability analysis, we find 


[2 


i sin 6 


0 - kAx 


( 2 . 10 ) 


( 2 . 11 ) 


( 2 . 12 ) 


Again, for stability, Reo < 0, and we now require that A < 0. 

In summary, for one-sided spatial difference approximations 
we have the following result: If 3u/3x is approximated by the 

backward-difference operator (2.6), then the resulting ordinary 
differential equation (2.7) will be stable if and only if A > 0 
(i.e., the wave travels to the right). Conversely, if 3u/3x is 
approximated by the forward-difference operator (2.11), the result- 
ing ordinary differential equation (2.7) will be stable if and 
only if A < 0 (i.e., the wave travels to the left). In general, 
no conventional backward, forward, or unsymmetric operator such as 


“2uj . - 3u 4 + 6u, - u, 

tl i±i it! + 0 (Ax 3 ) 

6Ax 


(2.13) 


8 



will yield an ordinary differencial equation which is simultaneously 


stable for both positive and negative eigenvalues . Although this 
statement is justified in Appendix A, its correctness is apparent 
from the fact that any noncentered spatial-difference operator will 
yield an eigenvalue (see, e.g., Eq. (2.9b)) with a nonzero real 
part whose coefficient is the eigenvalue X. Hence, the real part 
cannot satisfy Reo < 0 for both positive and negative X. 

Returning to the system (2.1), it is cltoar that if a single 
noncentered difference operator is used to approximate 3F/3x 
when the eigenvalues of the Jacobian matrix A are of mixed sign, 
then the resulting time-continuous method will always produce a 
numerical instability. 

3. ONE-DIMENSIONAL EQUATIONS OF GASDYNAMICS 

In one spatial dimension the inviscld equations of gasdynamics 
can be written in the conservation-law form (2.1) where 



V 


" m 

u - 

m 

F(U) - 

(m 2 /p) +p 




_(e + p)m/p_ 


and where m - pu. The primitive variables of (3.1) are the den- 
sity p, the velocity u, and the pressure p. The total energy 
per unit volume, e, is related to the internal energy per unit 
mass, e, by 

e ■ pc + pu 2 /2 ■ pc + m 2 /(2p) . (3.2) 
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The system is completed with an equation of state 

p - p(p,t) . (3.3) 

For the case of a perfect gas, 

p ■ (y - l)pe . (3.4) 

which can be rewritten using (3.2) as 

P - (Y - l)[e - m 2 /(2p)j , (3.5) 

where y is the ratio of specific heats. 

By using (3.5), the flux vector F(U) can be rewritten as 
" m 

F(U) » (y - l)e + (3 - y)m 2 /(2p) . (3.6) 

_yem/p - (y - l)m 3 /(2 p 2 ) __ 

The Jacobian matrix A * 3F/3U is easily computed and found to be 
0 1 0 
(y - 3'u 2 /2 (3 -y)u Y - 1 • (3.7) 

[(Y “ l)u 3 - yeu/p ye/p - 3 (y ~ l)u 2 /2 yu _ 

The eigenvalues of A are 

Xj«u, Xj • u - c , (3.8) 

where c ■ (Yp /p)^ 2 is the local speed of sound. For subsonic 
flow |u| < c, and the eigenvalues are of mixed sign since u + c 
and u - c are of opposite sign. 

The inviscid equations of gasdynamics have the rather remark- 
able property that if the equation of state has the functional 
form 

p - pf(c) , (3.9) 
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then the nonlinear flux vector F(U) is a homogeneous function of 
degree one in U; that is F(oU) ■ oF(U) for any value a. The 
equation of state (3.4) is clearly a special case of (3.9) and the 
fact that F(U) is a homogeneous function of degree one is obvious 
by inspection of the flux vector (3.6). By application of Euler’s 
theorem on homogeneous functions (see, e.g., [5]) there follows 

F - AU , (3.10) 

where A is the Jacobian matrix 3F/3U. One can readily verify 
the above equality by using Eqs. (3.7) and (3.1a) and making the 
indicated matrix-vector multiply. The flux vectors in two and 
three spatial dimensions also have the homogeneous property. 

If F satisfies the homogeneous property and A has a com- 
plete set of linearly Independent eigenvectors, then the flux vec- 
tor F can be split into subvectors, each one of which is asso- 
ciated with a tailored set of eigenvalues. In particular, the 
eigenvalues associated vith one subvector can be all positive, 
those associated with the other all negative. These subvectors can 
then be differenced Individually with an appropriate one-sided 
scheme in conservation-law form. The details for the one-dimensional 
case are outlined in the following section. 
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4. FLUX VECTOR SPLITTING FOR THE ONE-DIMENSIONAL EQUATIONS 

OF GASDYNAMICS 

Consider Eq. (2.1) with U and F defined by Eq. (3.1). The 
flux vector F(U) has the homogeneous property defined in the pre- 
ceding section and consequently F can be split into two pnrts 
as [2], [3] 

F = F* + F - , (4.1) 

where F + corresponds to the subvector associated with the posi- 
tive eigenvalues of A, and F~ corresponds to the negative eigen- 
values. This splitting is derived as follows. By virtue of (3.10) 
and (2.3), 

F = AU = QAQ _1 U , (4.2) 

where the diagonal elements of A are given by (3.8). 

Any eigenvalue X ^ can be expressed as 

A* - 4 + *£ > (4.3) 

where 

2 ’ “ 2 » (4 ‘ 4) 

so that if Aj, > 0 then A"£ - A^, ■ 0, with the converse result 

for Aj£ < 0. 

Using the above formulas, we split the diagonal matrix 

A " A + + A~ , (4.5) 

where A + and A” have as diagonal elements and A^, respec- 
tively. Equation (4.2) can be rewritten as 
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F = Q(A + + A")Q _1 U 
= (A + + A")U 


** 

where 

= F + + F~ , 

(4.6) 

m* 


A + = QA + Q~ 1 , A~ = QA'Q" 1 

(4.7) 



F + = A + U , F* = A'U , 

(4.8) 


and 

A = A + + A" . 

(4.9) 


The eigenvalues of A + are nonnegative and those of A" are non- 
positive. For the inviscid gasdynamic equations, the matrices Q 
and Q” 1 are given by 


Q - MT , Q" 1 * T -1 ^ 1 , (4.10) 

where M and T and their inverses are given in [2] for one and 
two space dimensions and in [6] for three space dimensions. 

The eigenvalues given by Eqs. (3.8) are split according to 
Eqs. (4.3) and (4.4) into 



The corresponding subvectors F + and F~ for the special case 
0 s u £ c are 
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F » 


_P_ 

2Y 


F s -K- 

2y 


2yu + c - u 

2(y-l)u 2 + (u+c) 2 

(Y - Du 3 + (u - ^ ) 3 + II- T . ) (u + c)c ; 

Ky } 2 2(y - 1) 


u - c 
(u - c) 2 

(u- c) 3 (3- V) (u- c)c 2 

rt ’ n / -« \ 


2(y - 1) 


or, if u > c, 


(4.12a) 


F+ = F , F" = 0 , 


(4.12b) 


where F is (3.6). The subvectors (4.12a) can be obtained, by a 
tedious calculation, directly from (4.6) or from the generalized 
flux vector (4.19) given at the end of this section. 

The eigenvalue splitting (4.4) is not unique and other split- 
tings into positive and negative parts are possible. For example, 
consider the splitting given by ( 3 ] : 


,+ _ u + u | . - _ u - 

A 1 " 2 A 1 2 

*2 " + C X 2 = A T 



(4.13) 


which satisfies (4.3) and the sum + = 1$, gives the physi- 

cal eigenvalues. 

A more general class of splitting is given by 

\ = X l + + • ' ‘ ’ (4 ’ 14) 
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with corresponding subvectors 


F - F a + F B + . . . , 


(4.15) 


where X^ is not necessarily split into positive and negative 
parts. For example, the sound wave or pressure term contribution 
to the eigenvalue could be split from the flow velocity. In the 
general notation (4.15) with a = u and $ = c, one has 

Xj = u Xj = 0 , | 


,u 

X 2 = u 


(4.16) 


% U * Vw- 

= U A. = “C , 


and the subvectors are given by 


F = uU = u I m 


F - P 


(4.17) 


This splitting, in two and three dimensions, has been used in 
parabolized Navier-Stokes calculations (7]. 

Since several splittings are possible, it is convenient to 
define a "generalized" flux vector from which any split subvector 
can readily be computed. The generalized flux vector is defined by 


JT = QAQ-^U , 


(4.18) 


where A is the diagonal matrix 





whose eigenvalues A^ are arbitrary. A direct calculation for 
the one-dimensional gasdynamic equations yields 

2 (y - 1) A j + A 2 + A 3 
2(y-l)AiU + A 2 (u + c) + A 3 (u-c) 




(y- DAjU 2 + — (u+ c) 2 + (u- c) 2 + w 


(4.19) 


where 


(3 - Y)(A 2 + A 3 )c J 


w 


(4.20) 


2(y - 1) 

The vector has a rather striking structure. One can 

A 

easily verify that if the A^ are replaced by the physical eigen- 
values (3.8), then (4.19) reduces to the physical flux vector (3.6) 
which, of course, must follow by Eq. (4.2). Flux formulas for any 
splitting follow directly by inserting the appropriate split 
eigenvalues (4.14) into Eq. (4.19). In particular, to arrive at 
the splitting defined by Eq. (4.11), F + follows directly by 
inserting A+ into Eq. (4.19), and F" follows by inserting A~ 
into Eq. (4.19). The matrices A + and A” can be obtained from 
Eq. (4.7). The other splittings are obtained in a similar way. 

The generalized flux vectors ^11 • *^11 for two and three 
spatial dimensions are given in Appendix B. From these generalized 
vectors, the flux vectors for any eigenvalue splitting can easily 
be computed for the inviscid equations of gasdynamics. 

In concluding this section we remark that, in general, 

A + 4 3F + /3U and A“ 4 3F“/3U. However, in all the numerical tests 
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that we have made, 9F + /9U does have positive eigenvalues and 
9F~/9U has negative eigenvalues; these roots, however, are not 
Identical to those of A + and A ~. 

5. ALGORITHMS FOR ONE SPACE DIMENSION 


In this section we illustrate several numerical algorithms 
that can be constructed for the one-dimensional equations of gas- 
dynamics by use of flux vector splitting. 


Explicit Methods 

MacCormack's scheme [8] for the one-dimensional system of 
conservation laws (2.1) is 


..n+1 n 

u j " u j - At • 

A 

U n+1 - 1 (u n+l + _ At jU_ 

J 2 V 2 Ax 


(5.1) 


(5.2) 


where Uj denotes the finite-difference approximation to U, 
f” * F(u”), etc., and the forward and backward difference oper- 
ators are defined by Eqs. (2.11) and (2.6). 

Since the predictor (5.1) is one-sided (upwind), the corrector 
(5.2) can be modified as 


V 2 F? 

n+1 1 „ n „n+lv At * j 

u, ■ 2 (u j + > - t sr - 1 


V F 
At x 


n+1 


(5.3) 


j 2 j j 2 Ax 2 Ax 

to obtain a completely upwind second-order scheme [1]. A necessary 
local condition for the stability of the scheme Eqs. (5.1) and (5.3) 
is that all the eigenvalues of the Jacobian matrix A be positive. 
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MacCormack’s scheme can be modified by using a forward differ- 
ence in the predictor, and a backward difference in the corrector. 
Likewise, the upwind scheme (5.1), (5.3) can be altered by replac- 


ing V by A and V 2 by -A 2 . In this case a necessary condition 
for the stability of the altered scheme is that all the eigenvalues 
of the Jacobian matrix be negative. 

The eigenvalue splitting (4.11) or (4.13) of the previous 
section can be used so that a split upwind version of MacCormack’s 
scheme can be used when the eigenvalues are of mixed sign, that is, 


in a subsonic region. The split upwind algorithm is 

-T7 _ A X< F i> n 

• n+1 - U? - At - At , 

j Ax Ax 

r.+\n+] 


U 


i 


(5.4) 


U n+1 . 1 (u n + u n +i } _ At / 7 x (F J>_ VV 

j 2 ' U j U j J 2 \ Ax Ax 

, at MflO 

2 y Ax " Ax f ' 


) 


(5.5) 


According to linear stability theory, the scheme (5.4) and (5.5) 
is stable if and only if 


|**|At/Ax < 2 

for all eigenvalues A“ where A^ » A+ + A^ are the eigenvalues 
of the Jacobian matrix. 

The MacCormack scheme, (5.1) and (5.2), is a symmetric scheme 
in the sense that the grid point cluster is symmetric about the 
center point j at the completion of the corrector step. A sym- 
metric (explicit) second-order scheme has a predominantly lagging 
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phase error and an upwind scheme has a leading phase error [1]. 
The opposite phase error of the symmetric and the upwind schemes 
suggests that a considerable reduction of phase error would occur 
if the two schemes were alternated on successive time steps. A 
temporal switching of schemes is the basis of Fromm's method of 
zero-average phase error [9]. 


Implicit Methods 

A noniterative implicit finite-difference scheme, for a one- 
dimensional system of conservation laws is [2] 


( T , 6At , . n\ n At 5 

1 + 1 + £ 6 x A jj AU j ~ 1 + C 6 x F j + 1 + 


Ml 0-1 

t 4U j > 


(5.6) 


where I is the identity matrix, AU n * U n+1 - U n , A is the 
Jacobian matrix, and <5 X is an appropriate spatial difference 
operator.* In general, the spatial derivative approximations on 
the left- and right-hand sides of (5.6) can be different. The 
parameters 6,5 determine the particular time-differencing 
approximation used. Scheme (5.6) includes three well-known 


implicit formulas 


e 


l , 


5 » 0 trapezoidal formula; 

5 * 0 backward Euler; 

1 


0*1, 5 * y three-point backward. 

For a more general formulation that includes all linear multistep 
(time-differencing) methods see [10]. 


*In Eq. (5.6) and in similar equations throughout this paper. 


notation of the form (I + Atfi A)AU denotes AU + At<5 (AAD). 

X X 



With use of the split flux vectors (4.7), one-sided spatial 
difference approximations are possible. For example, 

t 1 + TTJ (vjl n + Vjr)] 4U ” 

■ + s x r j |n ) + rfr-r 1 • 


(5.7) 


where 


and 


6 b F 
x j 


<S f F 
x j 


3Fj - 4Fj., + Fj. 2 
2Ax 


-3 Fj + 4F j+1 - F j+2 
2Ax 


(5.8) 


(5.9) 


are second-order accurate one-sided difference operators. 

The splitting F * F + + F allows an approximate factoriza- 
tion of the left-hand side of (5.7) into the product of two 
operators as 


(5.10) 


This scheme is implemented by the sequence 

( J + rrr Vjl"K ' irh-W* ? T * ‘ft 1 ") + rh AV T' ■ 

(5.11a) 


( I + rrr *x^l n ) 40 j ■ in * • 


1 • ,n a 4 < ,n 

u j - u j + Au j • 


(5.11b) 

(5.11c) 
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Unlike Eq. (5.7), the solution of Eq. (5.10) does not require 
the solution of a block tridiagoual system because both Eqs. (5.11a) 


and (5.11b) lead to block bidiagonal systems. For example, by 
writing out (5.1ib) 

t 1 - rf? zi (a ?°K ' 10 J • rh s - < 5 - 12 > 


v.e see that the solution is achieved by a right to left sweep 
(decreasing j) with the inversion of the 3 * 3 matrix in the 
parenthesis on the left required at each mesh point. The eigen- 
values of this matrix are greater than or equal to unity and con- 
sequently the matrix is nonsingular for any An. 

The scheme (5.10) is second-order accurate, dissipative, and 
unconditionally stable for 6*1, C * 1/2 (according to linear 
theory). In one spatial dimension, computational efficiency can be 
lost in comparison to Eq. (5.6) with <5 X a three-point central 
operator. This is chiefly because A , A , F , and F are costly 
to form. In multidimensions, however, an advantage is achieved by 
avoiding the solution of block tridiagonal systems. 


6. NUMERICAL EXPERIMENTS IN ONE DIMENSION 


The numerical solution of a one-dimensional shock-tube flow 
was chosen to judge the viability of the numerical algorithms given 
in the previous section. As a model problem, consider a tube of 
large extent in which a diaphragm separates a perfect gas at rest 
with different static pressures but at a uniform temperature. 
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With rupture of the diaphragm, an expansion propagates into the 
high-pressure gas, while a shock wave, followed by a contact sur- 
face, propagates into the low-pressure gas. Details of this flow 
are described in standard texts (e.g., Liepmann and Roshko [11]). 

In our calculations, the initial pressure ratio across the 
diaphragm is taken as 5 to 1. The solution results for various 
methods are shown in Figs. 1 to 5 in terms of the nondimens ional 
density, p/p Q , where p 0 is the initial high-density gas. In all 
cases, the same spatial grid and time step are used and 
At/Ax - 0.2. 

The results for the explicit upwind scheme (5.4), (5.5) are 
shown in Fig. 1. Also shown are the exact locations of the shock 
and contact waves and exact constant density levels. The Initial 
location of the diaphragm is taken at x ■ 3.0. Overall, the 
numerical accuracy is good and although the contact wave is 
smoothed out, the overshoots are moderate. For comparison pur- 
poses, the results for this flow obtained using the conventional 
MacCormack scheme are shown in Fig. 2. The accuracy of the two 
schemes is comparable, with the exception of a large spike in 
density as a result of start up. MacCormack has shown that the 
addition of a dissipation term, especially in expansion regions, 
can control such spikes [12]. We did not program this version, 
however, in order to illustrate the effectiveness of alternating 
the centered and upwind schemes. 
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In Fig. 3 we show Che solution obtained using the upwind scheme 
(5.4), (5.5) to advance the solution for odd values of the time 

index n and using the MacCormack scheme (5.1), (5.2) for the even 

values of n. The results of this combined algorithm are clearly 

superior to the application of either of its constituents. The 

overshoots are much reduced, and the jumps are crisper. 

Results for the implicit upwind scheme (5.10), are shown in 
Fig. 4. Here, trapezoidal time'differenclng was used. Again the 
results are good and quite comparable to those obtained with the 
explicit procedure. Finally, in Fig. 5 we show the results obtained 
from the "conventional" implicit algorithm using centered spatial 
differencing and three-point backward time-differencing. A small 
amount of fourth-order numerical dissipation was also added [13]. 
Overall, the results are again quite comparable. 

7. FLUX VECTOR SPLITTING IN TWO SPACE DIMENSIONS 


In two spatial dimensions, a hyperbolic system of conservation 
laws has the form 


where for the 


U - 


Jt Jx iy 0 * 


invlscid gasdynamics equations 


m ■ 


m « 


* « 

P 


pu 


pv 

m 

, F - 

pu 2 + p 

» G - 

puv 

n 


puv 


pv 2 + p 

e 


u(e + p) 


v(e + p) 

at m 


m m 




(7.1) 


(7.2) 


■i 

■j 
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where m = pu, and n » pv. The primitive variables of (7.2) are 
density p, velocity components u and v, pressure p, and total 
energy per unit volume e. The equation of state is 


P ■ (Y - 1)| 

This system can be rewritten in quasi-llnear form as 


[•-$*£)] 


8U 4. A 3U 4. * 3U 


0 , 


(7.3) 


(7.4) 


where A and B are the Jacobian matrices 

A - 3F/3U , B - 3G/3U , (7.5a,b) 

As in the previous section, we use the fact that F(U) and G(U) 
are homogeneous functions of degree one in U and consequently 

F - AU , G - BU . (7.6a„b) 

For the inviscid equations of gasdynamics, the matrices 
A and B can be diagonalized as 


Q _1 AQ 


O 


u + c 


u - c. 


(k, ■ 1, k, 


0 ) 


(7.7) 


Q-lBQ 


O 


v+ c 


v - c. 


(k. ■ 0, k, » 1) 


(7.8) 


where c is the local speed of sound. The matrices Q and Q" 1 , 
as defined in Appendix B, are functions of the two parameters 
kj,k 2 and of the dependent variables. The values of (kj,k 2 ) 
indicated in the parentheses of (7.7) and (7.8) are the valuer 
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Chat diagonalize A and B. Since A and B dc not commute, they 
cannot be diagonalized by the same similarity transformation. 

The flux vectors F(U) and G(U) for the inviscid equations oi. 
gasdynamics can be split into subvectors, each of which depends on 
eigenvalues of the same sign exactly as in the one-dimensional 
case. A generalized flux vector 

& u - QAQ-ty , (7.9) 

analogous to (4.19) tor one space dimension, is given in Appendix B 

for two space dimensions. (For completeness, a three-dimensional 

version is also given.) By using the generalized flux vector, one 
+ + 

can compute F“, G~ for any desired eigenvalue splitting. For 
example, 

F^ - ^(k, - 1, k 2 - 0, x{, Xj, X+) (7.10) 

where as given by Eq. (B9' of Appendix B is evaluated using 

the particular values of the parameters (kj,k 2 ) and (ij, X 3 , X u ) 
indicated in the parentheses of Eq. (7.10). The positive eigen- 
values X+, A*, X+ are defined by Eq. (4.4), where Aj - u, 

X 3 ■ u + c, X 4 • u - c. 

8. ALGORITHMS FOR TWO SPACE DIMENSIONS 

Just as in the one-dimensional case, the split subvecto.* 
forms allow construction of novel numerical difference schemes. 
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Explicit Methods 

The natural extension of the upwind scheme (5.4), (5.5) to two 
spatial dimensions is 

v T.l - u ”.k - "*(VV> n + V*},/] - v y[ a y (G j,k ) ” + 7 y^ G j,k> n ] 

( 8 . 1 ) 

<1 ■ 1 {c + u ".k - '-x[‘x<^,k> n+ ‘ + 7 x< F I,k> n+1 ] 

- ''y[ a y <G j,k >n+l + V G l,k>” + '] * ‘ 4 x (r j,k> n ] 


'[ v y <G j,k )n - 4 y <G j,k >n ]} ' 


( 8 . 2 ) 


where x * jAx, y * kAy, and v x = At/Ax, * At/Ay. Although 
this upwind version of MacCormack’s scheme requires considerably 
more work than the conventional scheme, it is a more robust algo- 
rithm if the solution exhibits large spatial gradients. In addi- 
tion, as in the one-dimensional case, a very effective algorithm 
is obtained when the upwind scheme is alternated with the conven- 
tional MacCormack scheme on successive time steps. 

If only first-order time accuracy is required, a simple 
explicit scheme is given by 


n+i ,.n 
J j,k u j,k 


Ilf' ■ uj v - 4t(^F + + jV + JyG + + 6*G~) |”_ k , (8.3) 


where 5^ and 5^ are defined by (5.8) and (5.9). Such a scheme 
is practical for steady state problems and it could be the basis 


of a point relaxation algorithm. 
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Implicit Methods 

The natural extension of (5.7) to two spatial dimensions is 

fl + -r—y (v A + | n + A AT , T + V B + | n + A B~ | "SI AU 1 ? 

L 1 + 5 V X j ,k 1 x j , k 1 y j.k 1 y j.k 1 f J j,k 

= -( . iv )(sV , I n + 6 f fT . | n +6 b G'!‘ , | n 6 f gT v | n ^ + (rlrV U "”i • 

Vl + C/Vxj.k 1 x j ,k 1 y j.k 1 yj.k 1 ) \L+K) j.k 

(8.4) 

The left-hand side of the scheme (8.4) can be factored into 
the product of two operators as 

i 1 + h (’x A !, k i i 1 ” • + VU 1 ' ")] i 1 + h ( 4 * A J ,k ' ° ' + Vi , fc I ”)] SU j ' RHS <8 ' 4> ' 


1 + 5 

and "RHS (8.4)" denotes the right-hand side of (8.4). This 
scheme can be implemented by the sequence 


(8.5) 



lu J,k - 

RHS (8.4) , 

(8.6a) 

1 

t 1 + h (Vj,k |n + VIA)} 

II 

C f-j 
<] 

4U U • 

(8.6b) 

> 


u n+1 = 
j.k 

“j.k + ‘“j.k' 

(8.6c) 

; 


Compared to the class of centrally-differenced implicit 
schemes [2, 14, 15] the algorithm (8.6) has both advantages and 
disadvantages. Equation (8.6a) requires the solution of a sparse 
block lower diagonal matrix and Eq. (8.6b) requires the solution 
of a sparse block upper bidiagonal matrix. Consequently, the com- 
putational inversion work is much less than that of solving two 
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block tridiagonal matrix sequences, as a conventional central- 
differenced algorithm would require. Moreover, in three dimensions 
the plus-minus split subvectors can still be approximately factored 
into just two factors — a sparse upper block-triangular matrix and a 
sparse lower block-triangular matrix. In three dimensions, the use 
of central spatial differences requires the inversion of three 
block tridiagonal sequences. The upwind differences are also 
dissipative, so it is not necessary to add higher-order dissipa- 
tion terms. 

On the other hand, twice as many Jacobian matrices and flux 
vectors have to be formed with the plus-minus splitting. Further- 
more, these are more involved to form than the usual Jacobians, 

although with careful programming certain terms in A“ and B ± as 
+ + 

well as F~ and G can be formed simultaneously. 

Other difference schemes can be formulated that use the plus- 
minus flux vector splitting. An implicit second-order accurate 
scheme is given by 

[ I + h (vI,kl" + ' Vi.k |n+ vl.kl")] [ I+h Vj,k!"K,k- RHS «•*> • 

( 8 . 7 ) 

This factorization does require a block tridiagonal inversion in 
the y-direction. Consequently, viscous terms in the y-direction 
can readily be included into the implicit operator. 

A semi- implicit, first-order accurate scheme in time is 
obtained from Eq. ( 8 . 7 ) by dropping the factor (i + hA x Aj ^1°) and 
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letting 0=1, £ = 0. Indeed, once the flux vectors are broken 
into subvectors, a large number of difference schemes can be devised 
to achieve possible advantages in numerical accuracy, robustness, 
computational efficiency, and storage. 

9. NUMERICAL EXPERIMENTS IN TWO DIMENSIONS 

Numerical calculations of model problems in two dimensions 
have been used to verify the stability and practicality of the 
algorithms of the last section. For example, the explicit algo- 
rithm (8.1) was tested on a square uniform grid with periodic bound- 
ary conditions. Waves were followed in time for arbitrary (non- 
physical) initial data. Although no results are shown, it was 
noted that the upwind scheme required about 3 times more computa- 
tional time than standard MacCormack scheme. 

The implicit algorithm (8.5), was tested on a simple biconvex 
airfoil with linearized boundary conditions. A typical transonic 
airfoil solution result is shown in Fig. 6 for a free-stream Mach 
number of 0.84 and a body thickness ratio of 11.4. The eigenvalue 
splitting (4.4) was used. Also shown are the calculated results 
obtained from the conventional implicit algorithm using central 
differencing. Both calculations use the same grid and boundary 
conditions. The grid is clustered in x and y and uses 50 * 28 
points. The results shown are in good agreement for this coarse 
grid. 
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The numerical calculations illustrated two weaknesses of the j 

upwind scheme that have now been essentially corrected. Whenever J 

] . 

an eigenvalue changes sign, it is either suddenly set to zero or ; 

is suddenly nonzero. Elements of any one subvector are suddenly j 

changed and the local accuracy of the difference approximation can r 

suffer. In Fig. 6 one notices a small oscillation in the data at 

i : 

the sonic line where the (u - c) eigenvalue changes sign. This ' 

oscillation would actually appear much worse if it were not for the i 

! ; 

blending terms that are added to the eigenvalues to smooth out < ; 

sudden changes. An example of blended terms are the following: 



where are small positive numbers which smoothly approach zero 

as | | increases. We remark that nonconservative formulations 
are not afflicted with flux vectors that have discontinuous 
derivatives. 

As previously noted, A + j 3F + /9U, although the two matrices 
share eigenvalues of the same sign. For the shock-tube calcula- 
tions the time step was limited by accuracy considerations and no 
difficulty was encountered in using A + and A on the left-hand 

side of the implicit algorithm. In the steady state airfoil cal- 

+ + 

culation we did find that using A instead of 3F/3U Imposed 
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explicit-like time-step restrictions. Use of the true Jacobian 
resulted in a more robust algorithm. This is quite a different 
result than was found with convecticn-sound speed splittings [3] 
in which the similarity matrix had identical stability properties 
to the actual Jacobian. 

10. GENERAL CONSERVATION FORMS FOR HYPERBOLIC SYSTEMS 

The implicit algorithms developed in the previous sections 
were for a Cartesian coordinate system; however, computational 
fluid dynamical problems involve flows over (or through) arbitrarily 
shaped bodies. 

In this section we show that the previously derived algorithms 
can be made applicable to general flow fields. One method of 
handling complex geometries is to map the physical plane — for 
example, an airfoil in two spatial dimensions — into a rectangular 
computational plane. The desired transformation has the property 
that the airfoil surface is coincident with coordinate lines in the 
physical plane, and the airfoil surface lies along the boundary of 
the rectangular computational plane. The transform should also 
cluster grid points to regions where large spatial gradients occur. 
Since the actual numerical computation is carried out in a trans- 
formed rectangular plane with a uniform mesh, we review in this 
section the form of the transformed conservation-law equations and 
the corresponding difference algorithm. 
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It can be shown [16] that the conservation-law form (7.1) is 
retained under an arbitrary time-dependent coordinate transformation 
£ = £(x,y,t) , n - n(x,y,t) , t = t . (10.1) 

In particular, one obtains 


3U ,3F . 3G 
— + rrz + T~ ° 0 , 
3t 3£ 3n 


( 10 . 2 ) 

where U ■ U/J and the flux vectors F and G are linear combina- 
tions of the vectors of (7.2): 

F * (£ t U + £ x r + € y O /J (10.3a) 

G « (nU + n F + n G)/J (10.3b) 

t x y 


and 


J - IILlH 1 

3(x,y) 


C x n y * C y n x 


(10.4) 


is the Jacobian of the transformation. It is important to note 
that Cartesian components of velocity and momentum are retained in 
(10.2). The equations in three spatial dimensions are straight- 
forward generalizations of the above equations. 

As in the previous sections, we use the fact that F(U) and 
G(U) are homogeneous functions of degree one in U. As a 
consequence 

F(U/J) - F(U) - AU , 


F (U) 
J 


j 


G(U/J) - G(U) ■ BU 


(10.5a) 

(10.5b) 


where 


A - 3F(U)/3U , B «= 3G(U)/3U . 
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(Note that A and 3 are the Jacobians for Cartesian coordinates 
defined by Eq. (7.5).) 

Hence, Eqs. (10.3) can be rewritten as 


F = (k 10 I + k u A + k 12 B)U , ' 
G = (k 2p I + k 21 A + k 22 B)U , ( 


( 10 . 6 ) 


where k lQ » , etc., are the scale factors. The generalized 

flux vector defined in the appendix can be used to calculate 

a + * + 

F", G for any arbitrary eigenvalue splitting. This is apparent 
from Eqs. (B7), (B6), and (Bl) of Appendix B because the coeffi- 
cients k} and k 2 of Eq. (Bl) are arbitrary real numbers, which 
for the present application are taken to be the scale factors 

E , E , E , etc. Finally, since the conservation-law form is 
t x y 

retained by (10.2), the numerical algorithms of Sec. 7 are 
directly applicable to general conservation forms. 


11. CONCLUDING REMARKS 


A hyperbolic system of conservation laws whose associated 
Jacobian matrices have positive and negative eigenvalues can only 
be spatially differenced as a system with centered operators. 
However, splitting the flux vectors into subvectors whose associated 
eigenvalues are of the same sign allows use of one-sided (upwind) 
operators. 

In this paper, we have made use of the fact that flux vectors 
of the inviscid gasdynamic equations are homogeneous functions of 


33 


degree one to construct flux vector splittings. As a consequence, 
new explicit and implicit dissipative difference methods are 

devised which are more robust and computationally efficient than 

conventional spatially centered schemes. Preliminary computational 

experiments show that the new methods are feasible, although 

clearly both additional analysis and numerical testing on "realistic" 

problems are required. 

APPENDIX A: INSTABILITY OF ONE-SIDED SCHEMES FOR HYPERBOLIC 

SYSTEMS WITH EIGENVALUES OF MIXED SIGN 


Here we examine the hyperbolic system of equations 


i£ + k SSL m o 

3t A 3x U » 


where A is a constant matrix with both positive and negative 
eigenvalues. In Sec. 2 it was shown that approximating 9 X U with 
either V X U or A X U must always lead to instability for the time- 
continuous system of equations. In this Appendix, we argue that 
all conventional backward, forward, or biased (e.g., Eq. (2.13)) 
finite-difference approximations to 9 X U will be unstable if A 
has both positive and negative eigenvalues. 

As shown in Sec. 2, Eq. (Al) can be transformed into an 
uncoupled system of scalar wave equations of the form 

|a + A|H.°, (A2> 
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du 

dt 


where X is either positive or negative. With Introduction of a 
spatial-difference approximation, Eq. (A2) is replaced by a system 

of ordinary differential equations 

Mu + ? , (A3) 

where u ■ (u, u 2 , u 3 , . . . , u j , . . . ) t , uj - u(JAx), ? contains, 
known boundary data (if any) and M is a constant coefficient 
matrix. For example, if 8 X is approximated by V x with given 
boundary data on the left-hand boundary, then 


M - - f- 
Ax 


’ 1 


* m 

u o 

- 1 1 o 
-11 


0 

0 

-1 1 l 

* x 

0 

• • 

* * " Ax 

• 

o • • 

« • 


* 

-1 1 

m m 


• 


(A4) 


If 3 x is approximated by Eq. (2.13), then 


M - - 


6Ax 


-3 6 -1 


X 

-2 -3 6-1 0 


0 

-2 -3 6 -1 

-i- 

’ r 6 Ax 

0 

• • • • 


0 

o * • • * 

• • • • 


* 

• 

• m 


• • m 


(A5) 


In this latter example the values of M would change to accommo- 
date an alternative choice of differencing at the right-hand 
boundary . 
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The solution of Eq. (A3) Is (cf. [17, 18]) 

u ■ u D + particular solution . 

Mt 

If M has an eigenvalue with a positive real part, then e •*> • 
as t -*• <*> and the solution is unbounded. However, any conven- 
tional backward, forward, or biased differenced scheme leads to a 
matrix M with a nonzero real trace. The trace contains X as a 
multiplier and the sum of the eigenvalues of M equals its trace. 

Consequently, for either the case X positive or X negative, 

Mt 

e -*■ » as t • since at least one eigenvalue will have a 

positive real part. Thus, conventional backward, forward, or 

biased difference schemes are unstable. 

The above argument, showing instability, is not altered by 

boundary conditions since they will only affect a few end-point 

elements of the matrix M. For a large matrix (refined spatial 

grid), these few elements cannot alter the sign of the trace. Note 

that the use of central spatial differencing, which can be stable, 

leads to a matrix M with zero trace. We remark also that a non- 

conventlonal upwind differencing such as 

3u| . r J-l ' u J-2 

9x|j Ax * 

also has a zero trace and thus must be proved to be unstable by 
another argument. 
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APPENDIX B: GENERALIZED FLUX VECTOR FOR TWO AND THREE SPATIAL I - 


DIMENSIONS 

Two-Dimensional Case 

The system of conservation lavs (7.1) can be rewritten in 
quasi-linear form (7.4). Define a matrix P as 

P - k t A + k 2 B , (Bl) 

where k 1 and k 2 are arbitrary real numbers. The system (7.4) is 
hyperbolic at the point (x, y, t, U) if there exists a similarity 
transformation such that 



where the eigenvalues X^ are real and the norms of Q and Q ” 1 
are uniformly bounded. 

The formulas in the remainder of this appendix pertain to the 
inviscid equations of gasdynamics for a perfect gas. The Jacobians 
A and B are 4x4 matrices and Q and Q ” 1 can be written as 

Q « MT , Q“* - T^M -1 , (B3) 

where M, T, and their inverses are given in the appendix of [2]. 

In general, the elements of Q and Q“* are functions of the param- 
eters (kj,k 2 ) and the dependent variables. For example, the ele- 
ments (Q) 33 and (Q — 1 ) 3 3 are 
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1 


, 33 - (v + <* 2 > , (Q~ 1 ) 3 3 " T~ 

/2pc 


/2c 


[ck, - (y - l)u] , 


where u and v aie the x- and y-velocity components, p is the 
density, c is the local speed of sound, y is the ratio of spe- 
cific heats, and 

kj » k j / (k 2 + k 2) 1 / 2 , ic 2 ■ k 2 / (k 2 + k^) 1 / 2 . (B4) 

The eigenvalues A^ of P are 

Aj * A, » k,u + k 2 v , A 3 « Aj + c(k^ + k 2 ) 1 / 2 , 

A 4 ■ Aj - c(k^ + k 2 ) 1 / 2 . 

Formula (B2) can be rewritten as 

P - kjA + k 2 B - QAQ ’ 1 . (B6) 


(B5) 


Hence, the matrix A or B, or any linear combination, can be 
recovered from (B6). 

As in the one-dimensional case (see Sec. 4), it is convenient 
to define a generalized flux vector by 

- QAQ -1 U , (B7) 

where now the eigenvalues Ag of the diagonal matrix A are taken 
to be arbitrary. Although the matrix Q~* is rather complex, the 
product Q -1 U is simply 


Q“*U 


y 

o 



(B8) 
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Completing the computation for the generalized flux vector we 
obtain 



2y 


2(y- DAj + A 3 4 

2(y - 1)XjU + :' 3 (u + ckj) + ^(u- ckj) 

2(y- 1 )A 1 v + A 3 (v + cic 2 ) 4 \(v- ck 2 ) 

* 

(Y“ 1)Xj(u 2 + v 2 ) + ~Y [(u+ckj) 2 4 (v+ck 2 ) 2 ] 

\ 

4 ~2 l(u- ckj) 2 4 (v- ck 2 ) 2 ] 4 Wjj 


(B9) 


where 


<3 - y)U 3 + \)c 2 
W II " 2(y - l) 

and kj, k, are defined by Eq. (BA). The conventional flux vector 
F(U) is obtained from (B9) if kj * 1, k 2 ■ 0, and the X ? 's as 
given by (B5) are inserted in (B9). Likewise, G(U) is recove' ed 
if kj ■ 0, k 2 ■ 1. Formula (B9) can be used to obtain any flux 
vector splitting as described in Secs. 7 and 10. 


Three-Dimensional Case 

In three spatial dimensions, a hyperbolic system of conserva- 
tion laws has the form 


!« + *im + igm + M2). . o 

at 3x ay a* 

This system can be rewritten in quasi-linear form as 


au . . au 
aF + A H 




B f 4 C |S - 0 

ay a* 




(BIO) 


(Bll) 
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where A, B, C are the Jacobian matrices 


A = 


3F 


B = 


3G 


C = 


3U ’ 3U ’ 

The generaiization of (Bl) is 


P = kjA + k 2 B + k 3 C . 


3H 

3U 


(B12) 


(B13) 


The eigenvalues of P are 


*1 = X 2 = = kjU + k 2 v + k 3 w , 1 

- Xj + c(k-k) l/2 , X 5 * Xj - c(k-k) 1 / 2 ,) 


(B14) 


where 


k*k * k 2 + k 2 + k 3 . 


The matrices M, T, and their inverses, which are needed to compute 
Q and Q” 1 as defined by Eq. (B3), are given in [6]. The gener- 
alized flux vector for three space dimensions is defined by 

= QAQ _1 U , (B15) 

where the eigenvalues X^ of the 5*5 diagonal matrix A are 
again arbitrary. For the purpose of calculating we assume 


that Xj = X 2 


X 3 because this is the case for the physical 


eigenvalues of the matrix P defined by Eq. (B14) and for the 
eigenvalue splittings of interest. The product Q~*U is 

r pkj(y - 1)1 


Q"‘U 


pk 2 (y - 1) 
pk 3 (y - 1) 
c //2 
c/fi 

4 
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(B16) 



and the generalized flux vector is 


^II 


2(y-1)Aj + A 4 + A s 


2(y - 1) AjU + A 4 (u + ckj) + A 6 (u- ckj) 


2(y - l) AjV + A 4 (v+ ck 2 ) + A 5 (v- ck 2 ) 


P 

2y 


2(y- 1)AjW + A 4 (w+ck 3 ) + Ag(w- ck 3 ) 


(y - 1) Aj (u 2 + v 2 + vr) + -y [(u + ck^)*- + (v+ck,,)*- + (w+ckj)*-] 

a 5 

+ ~y [(u- ckj ) 2 + (v- ck 2 ) 2 + (w- ck 3 ) 2 ] + W in + P 


(B17) 


where 


W 


(3 - Y)(A 4 + A 5 )c 2 
III = 2(y -1) 

P = 2p (y - DAjkjCkjW - k 3 v) 


u, v, and w are the x-, y-, and z-velocity components, and 
kj * kj/Ck^ + k| + k^) 1 / 2 , k 2 » k 2 /(k 2 + k 2 + k 2 ) 1 / 2 

k 3 - k 3 /(k 2 + k 2 + k 2 ) 1 / 2 . 
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MacCORMACK SCHEME 

SOLUTION AT TIME = 1.0 
INITIAL PRESSURE RATIO = 5.0 
(At/Ax) = 0.2 

O NUMERICAL 



x 


FIG. 2. Shock-tube solution obtained with explicit MacCormnck 
scheme. 
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ALTERNATING MacCORMACK AND 
EXPLICIT UPWIND SCHEMES 


SOLUTION AT TIME ■ 1.0 
INITIAL PRESSURE RATIO = 5.0 

(At/Ax) = 0.2 



x 


HG. 3. Shock-tube solution obtained by alternating explicit 
upwind and MacCormack schemes. 


47 



IMPLICIT UPWIND SCHEME 
SOLUTION AT TIME » 1.0 



FIG. 4. Shock-tube solution obtained i'roin implicit upwind 
scheme . 
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IMPLICIT CENTERED SCHEME 


SOLUTION AT TIME -1.0 
INITIAL PRESSURE RATIO * 5.0 
<At/Ax) = 0.2 

SMOOTHING COEFFICIENT OF 0.03 



x 


FIG. 5. Shock-tube .solution obtained using implicit algorithm 
with central spatial differencing. 
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O SPLIT UPWIND SCHEME 

CENTRAL DIFFERENCING ALGORITHM 

O 



FIG. 6. Steady .state solution for 11.4% thick parabolic arc 
airfoil, a 0.84. 
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